/*
 * ProximitySearch.h
 *
 *  Created on: Jan 15, 2011
 *      Author: LyonsDesktop
 */

#ifndef PROXIMITYSEARCH_H_
#define PROXIMITYSEARCH_H_

#include "DistanceFunction.h"
#include "../../Datasets/Dataset.h"
//#include "../../Extern.h"

namespace CLEVER {
namespace RegionDiscovery {
namespace DistanceFunctions {
//! Provides spatial subdivision with a kD-tree for numeric datasets.
class KDSubdiv
{
public:
	struct Node
	{
public:
		std::vector<int> indices;
		int axis;
		double distance;
		NumericDataset * d;

		void GetNearestNeighbour(DistanceFunction<double> * distfunc, bool SkipDiag,
		                         NumericDataset * comp, int iExample,
		                         int & iBest, double & DistBest)
		{
			if (children.size() > 0) {
				//children[0].GetNearestNeighbour(distfunc, SkipDiag, comp, iExample, iBest, DistBest);
				//children[1].GetNearestNeighbour(distfunc, SkipDiag, comp, iExample, iBest, DistBest);
				double v = comp->GetValue(iExample, axis);
				if (v < distance) {
					children[1].GetNearestNeighbour(distfunc, SkipDiag, comp, iExample, iBest, DistBest);
					if (fabs(v - distance) < DistBest) {
						children[0].GetNearestNeighbour(distfunc, SkipDiag, comp, iExample, iBest, DistBest);
					}
				} else   {
					children[0].GetNearestNeighbour(distfunc, SkipDiag, comp, iExample, iBest, DistBest);
					if (fabs(v - distance) < DistBest) {
						children[1].GetNearestNeighbour(distfunc, SkipDiag, comp, iExample, iBest, DistBest);
					}
				}
			} else   {
				int dsize = indices.size();
				//for (int i = 0; i < dsize; i++)
				//{
				//	int index = indices[i];
				//	if (SkipDiag && index == iExample) { continue; }

				//	double tdist = distfunc->Distance(d, index, comp, iExample);
				//	if (tdist < DistBest || iBest == -1)
				//	{
				//		iBest = index;
				//		DistBest = tdist;
				//	}
				//}
				//if (iBest != -1)
				//{
				//	double d0, d1;
				//	distfunc->DistanceSSE(d, iBest, iBest, comp, iExample, d0, d1, true);
				//	DistBest = d0;
				//}
				for (int i = 0; i < dsize; i += 2) {
					//double d0, d1;
					//distfunc->DistanceSSE(base, index, comp, iExample,
					//	(iExample + 1 < dsize) ? (iExample + 1) : (iExample),
					//	d0, d1);
					int index = indices[i];
					int index2 = (i + 1 < dsize) ? indices[i + 1] : index;
					double d0, d1;
					distfunc->DistanceSSE(d, index, index2, comp, iExample, d0, d1, false);
					//double d0 = distfunc->Distance(base, index, comp, iExample);
					//double d1 = distfunc->Distance(base, index2, comp, iExample);

					if ((d0 < DistBest || iBest == -1) && (!SkipDiag || index != iExample)) {
						iBest = index;
						DistBest = d0;
					}
					if ((d1 < DistBest || iBest == -1) && (!SkipDiag || index2 != iExample)) {
						iBest = index2;
						DistBest = d1;
					}
				}
				//if (iBest != -1)
				//{
				//	DistBest = distfunc->Distance(d, iBest, comp, iExample);
				//}
			}
		}


		void Split(NumericDataset * d, int numaxis, int axis, std::vector<int> & indices, int depth = 0)
		{
			this->d = d;
			this->indices = indices;
			if (indices.empty()) { return; }

			int attr = indices.size();
			double t = 0;
			double min_t = DBL_MAX;
			double max_t = DBL_MIN;
			for (int i = 0; i < attr; i++) {
				double v = d->GetValue(indices[i], axis);
				t += v;
				if (v > max_t) { max_t = v; }
				if (v < min_t) { min_t = v; }
			}
			//double median = t / (double)attr;
			double median = (min_t + max_t) / 2.0;

			this->axis = axis;
			this->distance = median;

			std::vector<int> front_indices, back_indices;

			for (int i = 0; i < attr; i++) {
				double v = d->GetValue(indices[i], axis);
				if (v < distance) { back_indices.push_back(indices[i]); } else { front_indices.push_back(indices[i]); }
			}

			if (front_indices.size() > 10 || back_indices.size() > 10) {
				if (front_indices.size() > 5 && back_indices.size() > 5 && depth < 5) {
					Node f, b;
					f.Split(d, numaxis, (axis + 1) % numaxis, front_indices, depth + 1);
					b.Split(d, numaxis, (axis + 1) % numaxis, back_indices, depth + 1);
					front_indices.clear();
					back_indices.clear();
					children.push_back(f);
					children.push_back(b);

					indices.clear();
				}
			}
		}


		std::vector<Node> children;
	};

	Node root;

	void Populate(NumericDataset * src)
	{
		root = Node();
		std::vector<int> indices;
		int size = src->Size();
		for (int i = 0; i < size; i++) {
			indices.push_back(i);
		}
		root.Split(src, src->GetNumAttributes(), 0, indices);
	}


};

//! Abstract base class to locate neighbours.
template<typename T>
class ProximitySearch
{
public:
	//! Calculates the nearest neighbour in the dataset comp to the example base[iExample]
	__inline__ virtual int NearestNeighbour(Dataset<T> * base, Dataset<T> * comp, int iExample) = 0;

	__inline__ virtual int NearestNeighbour(Dataset<T> * base, DoubleArrayDataset * comp, int iExample) = 0;

	__inline__ virtual int NearestNeighbourIndices(Dataset<T> * base, Dataset<T> * comp, int iExample, std::vector<int> & indices, int extra_index) = 0;

	virtual void BuildDistanceMatrix(NumericDataset * D) = 0;

	__inline__ virtual int NearestNeighbourMatrix(std::vector<int> & indices, int iExample) = 0;

	virtual void Release() { };
};

template<typename T>
class LinearProximitySearch : public ProximitySearch<T>
{
public:
	DistanceFunction<T> * distfunc;
	bool SkipDiag;
	bool UseSSE;

	//std::vector<std::vector<double>> distmat;
	double * distmat;
	int distmat_pitch;
	NumericDataset * Dmat;

	void BuildDistanceMatrix(NumericDataset * D)
	{
		Dmat = D;
		int dsize = D->Size();
		distmat_pitch = dsize;

		//if (distmat) { delete[] distmat; }
		//distmat.clear();
		distmat = new double[distmat_pitch * distmat_pitch];
		for (int i = 0; i < dsize * dsize; i++) { distmat[i] = 0; }
		//#ifndef FRAMEWORKPROFILER
		//	#pragma omp parallel for schedule(dynamic, 64)
#pragma omp parallel for
		//#endif
		for (int i = 0; i < dsize; i++) {
			//std::vector<double> row;
			for (int j = i + 1; j < dsize; j++) {
				distmat[(i * distmat_pitch) + j] =
				    distmat[(j * distmat_pitch) + i] = distfunc->Distance(D, i, D, j);
				//row.push_back(distfunc->Distance(D, i, D, j));
			}
			//distmat.push_back(row);
		}
	}


	inline int NearestNeighbourMatrix(std::vector<int> & indices, int iExample)
	{
		int iNN = 0;
		double mindist = DBL_MAX;
		int dsize = (int)indices.size();
		int * baseindex = &indices[0];

		for (int i = 0; i < dsize; i++) {
			int index = baseindex[i];

			if (SkipDiag && index == iExample) { continue; }

			double tdist = distmat[(iExample * distmat_pitch) + index];
			//double tdist = distmat[index][iExample];
			if (tdist < mindist) {
				iNN = i;
				mindist = tdist;
			}
		}
		return iNN;
	}


	LinearProximitySearch(DistanceFunction<T> * distfunc, bool SkipDiag, bool UseSSE = false)
		: distfunc(distfunc), SkipDiag(SkipDiag), UseSSE(UseSSE), Dmat(0), distmat(0)
	{

	}


	void Release()
	{
		delete[] distmat;
		distmat = 0;
	}


	~LinearProximitySearch()
	{
		Release();
	}


	//! See: ProximitySearch::NearestNeighbour
	int NearestNeighbour(Dataset<T> * base, Dataset<T> * comp, int iExample)
	{
		if (UseSSE) { return NearestNeighbourSSE(base, comp, iExample); }

		int iNN = 0;
		double mindist = DBL_MAX;
		int dsize = base->Size();

		for (int index = 0; index < dsize; index++) {
			if (SkipDiag && index == iExample) { continue; }

			double tdist = distfunc->Distance(base, index, comp, iExample);
			if (tdist < mindist) {
				iNN = index;
				mindist = tdist;
			}
		}
		return iNN;
	}


	//! See: ProximitySearch::NearestNeighbour
	int NearestNeighbour(Dataset<T> * base, DoubleArrayDataset * comp, int iExample)
	{
		if (UseSSE) { return NearestNeighbourSSE(base, comp, iExample); }

		int iNN = 0;
		double mindist = DBL_MAX;
		int dsize = base->Size();

		for (int index = 0; index < dsize; index++) {
			if (SkipDiag && index == iExample) { continue; }

			double tdist = distfunc->Distance(base, index, comp, iExample);
			if (tdist < mindist) {
				iNN = index;
				mindist = tdist;
			}
		}
		return iNN;
	}


	//! See: ProximitySearch::NearestNeighbourIndices
	//indices: adds representatives, these representative are not exist in representative set of the previous iteration
	int NearestNeighbourIndices(Dataset<T> * base, Dataset<T> * comp, int iExample, std::vector<int> & indices, int extra_index)
	{
		if (UseSSE) { return NearestNeighbourSSEIndices(base, comp, iExample, indices, extra_index); }

		int iNN = 0;
		double mindist = DBL_MAX;
		int dsize = indices.size();

		for (int i = 0; i <= dsize; i++) {
			int index = (i < dsize) ? indices[i] : extra_index;
			if (index == -1)
				int temp = 1 + index;
			if (SkipDiag && index == iExample) { continue; }

			double tdist = distfunc->Distance(base, index, comp, iExample);
			if (tdist < mindist) {
				iNN = index;
				mindist = tdist;
			}
		}
		return iNN;
	}


private:
	//! Calls DistanceSSE in the distance functions and passes a pair per call.
	//! See: ProximitySearch::NearestNeighbour
	inline int NearestNeighbourSSE(Dataset<T> * base, Dataset<T> * comp, int iExample)
	{
		int iNN = 0;
		double mindist = DBL_MAX;
		int dsize = base->Size();

		for (int index = 0; index < dsize; index += 2) {
			//double d0, d1;
			//distfunc->DistanceSSE(base, index, comp, iExample,
			//	(iExample + 1 < dsize) ? (iExample + 1) : (iExample),
			//	d0, d1);
			int index2 = (index + 1 < dsize) ? (index + 1) : index;
			double d0, d1;
			distfunc->DistanceSSE(base, index, index2, comp, iExample, d0, d1, true);
			//double d0 = distfunc->Distance(base, index, comp, iExample);
			//double d1 = distfunc->Distance(base, index2, comp, iExample);

			if (d0 < mindist && (!SkipDiag || index != iExample)) {
				iNN = index;
				mindist = d0;
			}
			if (d1 < mindist && (!SkipDiag || index2 != iExample)) {
				iNN = index2;
				mindist = d1;
			}
		}
		return iNN;
	}


	inline int NearestNeighbourSSE(Dataset<T> * base, DoubleArrayDataset * comp, int iExample)
	{
		int iNN = 0;
		double mindist = DBL_MAX;
		int dsize = base->Size();

		for (int index = 0; index < dsize; index += 2) {
			//double d0, d1;
			//distfunc->DistanceSSE(base, index, comp, iExample,
			//	(iExample + 1 < dsize) ? (iExample + 1) : (iExample),
			//	d0, d1);
			int index2 = (index + 1 < dsize) ? (index + 1) : index;
			double d0, d1;
			distfunc->DistanceSSE(base, index, index2, comp, iExample, d0, d1, true);
			//double d0 = distfunc->Distance(base, index, comp, iExample);
			//double d1 = distfunc->Distance(base, index2, comp, iExample);

			if (d0 < mindist && (!SkipDiag || index != iExample)) {
				iNN = index;
				mindist = d0;
			}
			if (d1 < mindist && (!SkipDiag || index2 != iExample)) {
				iNN = index2;
				mindist = d1;
			}
		}
		return iNN;
	}


	//! Calls DistanceSSE in the distance functions and passes a pair per call.
	//! See: ProximitySearch::NearestNeighbourIndices
	inline int NearestNeighbourSSEIndices(Dataset<T> * base, Dataset<T> * comp, int iExample, std::vector<int> & indices, int extra_index)
	{
		int iNN = 0;
		double mindist = DBL_MAX;
		int dsize = indices.size();

		for (int i = 0; i <= dsize; i += 2) {
			//double d0, d1;
			//distfunc->DistanceSSE(base, index, comp, iExample,
			//	(iExample + 1 < dsize) ? (iExample + 1) : (iExample),
			//	d0, d1);
			int index = (i < dsize) ? indices[i] : extra_index;
			int index2 = (i + 1 < dsize) ? indices[i + 1] : extra_index;
			double d0, d1;
			distfunc->DistanceSSE(base, index, index2, comp, iExample, d0, d1, true);
			//double d0 = distfunc->Distance(base, index, comp, iExample);
			//double d1 = distfunc->Distance(base, index2, comp, iExample);

			if (d0 < mindist && (!SkipDiag || index != iExample)) {
				iNN = index;
				mindist = d0;
			}
			if (d1 < mindist && (!SkipDiag || index2 != iExample)) {
				iNN = index2;
				mindist = d1;
			}
		}
		return iNN;
	}


};
}
}
}

#endif /* PROXIMITYSEARCH_H_ */
